import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from scipy.io import loadmat
import pandas as pd

path = '/home/wdd/Desktop/ML_CODE/NeualWork/ex3data1.mat'
data = loadmat(path)
# print(data['X'].shape)
# print(data['y'].shape)

# 绘制100张图片
sample_idx = np.random.choice(range(data['X'].shape[0]), 100)
sample_images = data['X'][sample_idx, :]
# fig, ax = plt.subplots(nrows=10, ncols=10)
# for r in range(10):
#     for c in range(10):
#         ax[r, c].mathow(np.array(sample_images[10 * r + c].reshape((20, 20))))
#         plt.xticks(np.array([]))
#         plt.yticks(np.array([]))

fig, ax_array = plt.subplots(nrows=10, ncols=10, sharey=True, sharex=True, figsize=(12, 12))
for r in range(10):
    for c in range(10):
        ax_array[r, c].matshow(np.array(sample_images[10 * r + c].reshape((20, 20))).T, cmap=matplotlib.cm.binary)
        plt.xticks(np.array([]))
        plt.yticks(np.array([]))


# plt.show()

def sigmoid(z):
    return 1 / (1 + np.exp(-z))


# 代价函数,具备正则化参数
def cost(theta, x, y, learningRate):
    theta = np.matrix(theta)
    x = np.matrix(x)
    y = np.matrix(y)
    first = np.multiply(-y, np.log(sigmoid(x * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(x * theta.T)))
    reg = (learningRate / (2 * len(x))) * np.sum(np.power(theta[:, 1:theta.shape[1]], 2))
    return np.sum(first - second) / len(x) + reg


# 向量化梯度
def gradient(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)

    parameters = int(theta.ravel().shape[1])
    error = sigmoid(X * theta.T) - y

    grad = ((X.T * error) / len(X)).T + ((learningRate / len(X)) * theta)

    # intercept gradient is not regularized
    grad[0, 0] = np.sum(np.multiply(error, X[:, 0])) / len(X)

    return np.array(grad).ravel()


from scipy.optimize import minimize


def one_vs_all(x, y, num_labels, learning_rate):
    rows = x.shape[0]
    params = x.shape[1]

    # k * (n + 1) array for the parameters of each of the k classifiers
    # 设置10个分类器
    all_theta = np.zeros((num_labels, params + 1))

    # insert a column of ones at the beginning for the intercept term
    X = np.insert(x, 0, values=np.ones(rows), axis=1)

    # labels are 1-indexed instead of 0-indexed
    for i in range(1, num_labels + 1):
        theta = np.zeros(params + 1)
        # 对每一个分类器件分别进行优化,此时当前类别为1,非当前类别为0
        y_i = np.array([1 if laber == i else 0 for laber in y])
        y_i = np.reshape(y_i, (rows, 1))

        # minimize the objective function
        # 分别优化每一个分类器
        fmin = minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient)
        all_theta[i - 1, :] = fmin.x
    return all_theta


all_theta = one_vs_all(data['X'], data['y'], 10, 1)


# print(all_theta)


def predict_all(X, all_theta):
    rows = X.shape[0]
    params = X.shape[1]
    num_labels = all_theta.shape[0]

    # same as before, insert ones to match the shape
    X = np.insert(X, 0, values=np.ones(rows), axis=1)

    # convert to matrices
    X = np.matrix(X)
    all_theta = np.matrix(all_theta)

    # compute the class probability for each class on each training instance
    h = sigmoid(X * all_theta.T)

    # create array of the index with the maximum probability
    h_argmax = np.argmax(h, axis=1)

    # because our array was zero-indexed we need to add one for the true label prediction
    h_argmax = h_argmax + 1

    return h_argmax


if __name__ == "__main__":
    y_pred = predict_all(data['X'], all_theta)
